Axisymmetric core collapse simulations using characteristic numerical relativity 
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We present results from axisymmetric stellar core collapse simulations in general relativity. Our 
hydrodynamics code has proved robust and accurate enough to allow for a detailed analysis of the 
global dynamics of the collapse. Contrary to traditional approaches based on the 3+1 formulation of 
the gravitational field equations, our framework uses a foliation based on a family of outgoing light 
cones, emanating from a regular center, and terminating at future null infinity. Such a coordinate 
system is well adapted to the study of interesting dynamical spacetimes in relativistic astrophysics 
such as stellar core collapse and neutron star formation. Perhaps most importantly this procedure 
allows for the unambiguous extraction of gravitational waves at future null infinity without any 
approximation, along with the commonly used quadrupole formalism for the gravitational wave 
extraction. Our results concerning the gravitational wave signals show noticeable disagreement 
when those are extracted by computing the Bondi news at future null infinity on the one hand and 
by using the quadrupole formula on the other hand. We have strong indication that for our setup 
the quadrupole formula on the null cone does not lead to physical gravitational wave signals. The 
Bondi gravitational wave signals extracted at infinity show typical oscillation frequencies of about 
0.5 kHz. 

PACS numbers: 04.25.Dm, 04.40.-b, 95.30.Lz, 04.40.Dg, 97.60.Lf 



I. INTRODUCTION 

Supernova core collapse marks the final stage of the 
stellar evolution of massive stars. If the core collapse 
and/or the supernova explosion are nonspherical, part of 
the liberated gravitational binding energy will be emit- 
ted in the form of gravitational waves. According to esti- 
mates from numerical simulations, the total energy emit- 
ted in gravitational waves in such events can be as high 
as IQ-^Mqc^ 0,11 nil. Nonsphericity can be caused 
by the effects of rotation, convection and anisotropic 
neutrino emission leading either to a large-scale devia- 
tion from spherical symmetry or to small-scale statistical 
mass-energy fluctuations (for a review, see e.g., 0)- Su- 
pernovae have always been considered among the most 
important sources of gravitational waves to be eventually 
detected by the current or next generations of gravita- 
tional wave laser interferometers. If detected, the grav- 
itational wave signal could be used to probe the models 
of core collapse supernovae and to study the formation 
of neutron stars. 

Earlier studies of axisymmetric supernova core collapse 
were performed using Newton's law of gravity 
More recently, effects of general relativity have been in- 
cluded under the sirnplifying assumption of a conformally 
flat spatial metric 0, ■ In all existing works grav- 
itational waves are not calculated instantly within the 
numerical simulation, but they are extracted a posteri- 
ori using the approximation of the quadrupole formal- 
ism which links gravitational waves to the change of the 
quadrupole moment of the simulated matter distribution. 

In this paper we present first results of a project aimed 



at studying the dynamics of stellar core collapse by means 
of numerical simulations in full general relativity. The 
trademark of our approach is the use of the so-called 
characteristic formulation of general relativity (see P| for 
a review) , in which spacetime is foliated with a family of 
outgoing light cones emanating from a regular center. 
Due to a suitable compactification of the global space- 
time future null infinity is part of our finite numerical 
grid where we can unambiguously extract gravitational 
waves. This remarkable feature is the main motivation 
behind our particular choice of slicing and coordinates, 
which clearly departs from earlier investigations. We 
note, however, that we have not modeled the detailed 
microphysics of core collapse supernovae, which is be- 
yond the scope of the present investigation. Instead, we 
only take into account the most important features for 
both, the gravitational field and the hydrodynamics, and 
those will be introduced in the upcoming section. 

Characteristic numerical relativity has traditionally fo- 
cused on vacuum spacetimes. In recent years the field has 
witnessed steady improvement, and robust and accurate 
three-dimensional codes are nowadays available, as that 
described in which has been applied to diverse stud- 
ies of black hole physics (e.g. J,!]). In black hole space- 
times, only the geometry outside a horizon is covered by 
the foliation of light cones. This is is not the case for 
neutron stars or gravitational collapse spacetimes which 
must include a regular origin. Up to now, characteris- 
tic vacuum codes with a regular center have only been 
studied in spherical symmetry and in axisymmetry [l^ . 

The inclusion of relativistic hydrodynamics into the 
characteristic approach along with the implementation of 
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high-resolution shock-capturing (HRSC) schemes in the 
solution procedure was first considered by Papadopou- 
los and Font [T^ IT3 . . First applications in spherical 
symmetry were presented, dealing with black hole ac- 
cretion [l6j and the interaction of relativistic stars with 
scalar fields as simple models of gravitational waves. 
Axisymmetric studies of the Einstein-perfect fluid sys- 
tem were first discussed in |T^ . In this reference we 
presented an axisymmetric, fully relativistic code which 
could maintain long-term stability of relativistic stars 
and which allowed us to perform mode-frequency com- 
putations and the gravitational wave extraction of per- 
turbed stellar configurations. The core collapse simula- 
tions presented in the current paper are based on this 
code, which is described in detail in Ref. 

The paper is organized as follows: In Sec. |n] we de- 
scribe the mathematical and numerical framework we use 
in the simulations. Section IlIII deals with presenting the 
initial data of the unstable equilibrium stellar configu- 
rations we evolve. Section |W1 is devoted to discuss the 
numerical simulations, with emphasis on the collapse dy- 
namics. In Sec. we analyze the corresponding gravita- 
tional wave signals. A summary and a discussion of our 
results are given in Sec. lVIl Finally, tests to calibrate the 
code in simulations of core collapse are collected in the 
Appendix. 



II. FRAMEWORK AND IMPLEMENTATION 

We only briefly repeat the basic properties of our ap- 
proach here. The interested reader is referred to our pre- 
vious work jl8] for more details concerning the math- 
ematical setup and the numerical implementation. As 
described in Ref. 0, we work with the coupled system 
of Einstein and relativistic perfect fluid equations 



Gab — I^Tab , 


(1) 


VaT'^' = , 


(2) 


VaJ° = 0, 


(3) 



where Vo, as usual, denotes the covariant derivative. The 
energy-momentum tensor for a perfect fluid Tab takes the 
form. 



Tab = phUaUb + pgab- 



(4) 



Here p denotes the rest mass density, h = l + e + ^ is the 
specific enthalpy, e is the specific internal energy, and p 
is the pressure of the fiuid. The four- vector u", the 4- 
velocity of the fluid, fulfills the normalization condition 
gabu°'u^ = — 1. The four-current J" is defined as J" = 
pu". Using geometrized units (c — G — 1) the coupling 
constant in the field equations is k = Stt. We further use 
units in which Mq — 1. Moreover, an equation of state 
(EoS) needs to be prescribed, p — p{p,e), as we discuss 
in Sec. Hm below. 

Our numerical implementation of the field equations of 
general relativity is based on a spherical null coordinate 



system {u,r,9,4>)- Here, u denotes a null coordinate la- 
beling outgoing light cones, r is a radial coordinate, and 9 
and (p are standard spherical coordinates. Assuming axi- 
symmetry, (f) is a Killing coordinate. In order to resolve 
the entire radial range from the origin of the coordinate 
system up to future null infinity, we define a new radial 
coordinate x e [0, 1]. The radial coordinate r is a func- 
tion of the coordinate x, which can be adapted to the 
particular simulation. In this work, except where other- 
wise stated, we use a grid function r{x) = 100tan(^a;), 
for which the limit x ~* 1 corresponds to r cx3. More- 
over, in order to eliminate singular terms at the poles 
(9 — 0, tt) we introduce the new coordinate y = — cos9. 



A. The characteristic Einstein equations 

The geometric framework relies on the Bondi (radia- 
tive) metric 0| 



V 



ds-" = - 



2Ur'^e^''dud9 + r'^{e^'^d9'^ + e^^T sin^ 9d<lP). (5) 



We substitute the metric variables V, U, 7) by the new 
set of metric variables S, U, 7), 
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U 

7 



V -r 

U 
sin 9 ' 
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8111^9' 



(6) 
(7) 

(8) 



in order to obtain regular expression in the Einstein equa- 
tions in particular at the polar axis. The origin of the 
coordinate system r = is chosen to lie on the axis 
of our axisymmetric stellar configurations, where we de- 
scribe boundaries and falloff conditions for the metric 
fields. The complete set of Einstein equations reduces to 
a wave equation for the quantity 7 (see Eq. ^) and a 
hierarchical set of hypersurface equations for the quanti- 
ties (/3, U, S) to be solved along the light rays u — const. 
The particular form of these equations is explicitly given 
in Ref. ,1^. 



B. The relativistic perfect fluid equations 

The axisymmetric general relativistic fluid equations 
on the light cone, Eqs. Q and are written as a 
first-order flux-conservative, hyperbolic system for the 
state- vector U = (C/", ?7^, Uy, W^) = (T™, T"^, T^, J"). 
Following our previous work we have not included 
the metric determinant in the definition of the state 



vector. Explicitly, in the coordinates {x' 



x\ 
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{u,x,y, 



we obtain 



doU' 
doU 

d,Uy+d,F'y 

doU-^ + djF^^ 
The flux vectors are defined as 
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Sy , 
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(9) 
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(14) 
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and the corresponding source terms read 



= -{h-ii^)) ,^T^' + g^' St + n{9'''),c, (17) 

1 



Sa = --phu,u4g''')^a+p-{HV^)),,T', , (18) 



= -(ln(^/^)), J\ 



(19) 



wherein a comma is used to denote a partial derivative. 

The fluid update from time m" to it"+^ at a given cell 
i,j is given by 



/2j) 



Ax 

^(Gij + i/2 - Gij_i/2) 

AuSij , 



(20) 



where the numerical fluxes, F and G, are evaluated at 
the cell interfaces according to a flux-formula, the one 
due to Harten, Lax and van Leer (HLL) in our case [20l| . 
The characteristic information of the Jacobian matrices 
associated with the hydrodynamical fluxes, which is used 
in this flux formula, was presented elsewhere Il6|| . We use 
the monotonized central difference slope limiter by van 
Leer j2l| for the reconstruction of the hydrodynamical 
quantities at the cell interfaces needed in the solution 
of the Riemann problems. This scheme is second order 
accurate in smooth monotonous parts of the flow and 
gives improved results compared to the MUSCL scheme 
applied in (for an independent comparison, see [2^). 



C. Equation of state 

We use a hybrid EoS which includes the effect of stiff- 
ening at nuclear densities and the effect of thermal heat- 
ing due to the appearance of shocks. Such an EoS was 
first considered by Janka et al. (23i] . and has been used 
for core collapse simulations both, using Newtonian grav- 
ity 0, 13 and in general relativity under the assumption 
of conformal flatness 0, Q, 111 • 

In our EoS the total pressure consists of a polytropic 
part, which takes into account the contribution from the 



degenerate electron g well as the nuclear forces (at 
high densities), and a thermal part due to the heating of 
the material by a shock, p = Pp + pth- More precisely, 
the polytropic part follows the relation 



Pp 



Kip^\ for p < p„. 



K2P 



for p> Pn, 



(21) 



where we assume a nuclear density pn — 2 x 10^** g cm 
For a degenerate relativistic electron gas F — Fi„i = 4/3 
and K = 4.8974894 x 10^^^ [cgs]. To model the physi- 
cal processes which lead to the onset of the collapse, we 
reduce the effective adiabatic index from F to Fi setting 
Ki = K at the initial slice. Moreover, to model the stiffen- 
ing of the EoS at nuclear densities, we assume F2 = 2.5. 
The value of the polytropic constant K2 follows from the 
requirement that the pressure is continuous at nuclear 
density. The thermodynamically consistent internal en- 
ergy distribution reads 



Ti-lf^ 

K2 J^2 



\ forp<pn, 
+ E, for p > Pn- 



(22) 



The requirement that Cp is continuous at nuclear density 
leads to 



E 



(r2-Fi)^i 

(F2 - l)(Fi - 1) 



(23) 



For the thermal contribution to the total pressure, we 
assume an ideal fluid EoS 



Pth = (Fth - l)/oeth, 



(24) 



with an adiabatic index Fth = f describing a mixture of 
relativistic and non-relativistic gases. The internal ther- 
mal energy eth is simply obtained from 

eth = e - Cp. (25) 
We can summarize the EoS in a single equation: 



P 



K 1 



F - 1 



P^ + (Fth - l)pe 



(Fth-1)(F-Fi) 

(F2 - l)(Fi - 1) "^^^ ^' 



(26) 



where F and k change discontinuously at nuclear density 
Pn from Fi to F2 and ki to K2. For the sound speed Cg, 
we obtain 

/ic^ = i(Fpp + FthPth). (27) 



D. Recovery of the primitive variables 

After the time update of the state-vector of hydrody- 
namical quantities, the primitive variables {p,u^ ,Uy,e) 
have to be recomputed. The relation between the two 
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sets of variables is not in closed algebraic form. Using 
the hybrid EoS, such recovery is performed as follows: 
With the definition S"^'' = g'^'^T'^^T''^, we obtain [H 



P 



P 



5™ = l^i--l-ej{^i- + l + ej ( J")2 + (28) 

where in our null coordinate system g"" = 0. Let 
F{p,e) = |. From Eqs. ^ and the definition 

of the specific enthalpy h we obtain the 3 equations for 
the 3 unknowns F, p and e 



F{p,e) = ^L+{l + e)\ (29) 
F(p,e) = {Tth~l)e + G{p), (30) 
p - H{\ + e + F{p,e)). (31) 



In these equations we made use of the abbreviations 



L = 



(J" 



I = 



r - 1 

(r,,-i)(r-ri)_r_i 



(r2-i)(ri-i) 



■KP„ 



(32) 

(33) 
(34) 



From Eqs. (|29|1 - H31|I we deduce a single implicit equation 
for the rest mass density p 

/imp(p) (y)'-2^(l + e)-L = 0, (35) 
where we consider the internal energy as function of p 



.^j-(f-(i + GW)), 



(36) 



We solve Eq. (|35|) for p with a Ncwton-Raphson method. 



III. INITIAL DATA 

In the final stage of the evolution of massive stars, the 
iron core in the stellar center has a central density of 
about pc — 10^^ g cm~'^ when it becomes dynamically 
unstable to collapse. As the pressure of the degenerate 
relativistic electrons is by far the most important con- 
tribution to the total pressure, the pressure in the iron 
core can be approximated by a F = | polytropic EoS. 
In order to obtain an initial model for the iron core, 
we solve the Tolman-Oppenheimer-Volkoff equation [1^ 
with the above central density, which corresponds to 
pc = 1.62 X 10~* in code units (c = G = M© = 1). 

To initiate the gravitational collapse we set the adia- 
batic index Fi in the hybrid EoS H26|l to a value of 1.30, 
which mimics the softening of the EoS due to capture 
of electrons and due to the endothermic photodisinte- 
gration of heavy nuclei. The chosen value is within the 



interval range analyzed in previous studies of rotational 
core collapse based on Newtonian physics and on the 
conformal fiat metric approximation of general relativ- 

Since rotation is not included in our current implemen- 
tation, the equilibrium initial models of the iron core are 
spherically symmetric. Furthermore, in the evolution of 
these data during the phases of collapse, bounce, and be- 
yond, spherical symmetry is conserved. Therefore, since 
we are mainly interested in simulating core collapse as 
a source of gravitational waves, we add non-radial per- 
turbations on top of the spherical data. Our analysis 
is thus restricted to collapse scenarios where the effects 
of rotation are unimportant and in which stellar evolu- 
tion has led to asymmetries in the iron core, e.g. due 
to convection . The strongest gravitational wave sig- 
nals are expected for perturbations of quadrupolar form. 
Hence, we further restrict our analysis to this case, vary- 
ing the form and amplitude of the perturbation in the 
initial data. We note that the evolution of such data, 
however, can produce an arbitrary type of perturbation 
within the class of the imposed symmetry. 

We have classified the different models as follows: In 
case 21 the spherical model is unperturbed; in case 05 we 
prescribe a perturbation of the rest mass density 



Sp = Aps sin 



y 



(37) 



where pa denotes the spherical density distribution. Fi- 
nally, in case £ we prescribe a perturbation of the merid- 
ional velocity component 



Uy^As\n\ — ]y. 



(38) 



In the above two equations A is a free parameter describ- 
ing the amplitude of the perturbation, and R denotes the 
radius of the iron core {R = 1.4 x 10^ km). We note in 
passing that in [31 we already used a perturbation of 
the form £ to study quadrupolar oscillations of relativis- 
tic stars. We have further classified models 05 and £ 
according to the amplitude A of the perturbation (e.g. 
case £01 would correspond to an amplitude A — 0.1). 



IV. CORE COLLAPSE DYNAMICS 

This section deals with the description of the global dy- 
namics of our core collapse simulations. Relevant tests of 
the code which assess its suitability for such simulations 
are collected in the Appendix. 



A. Collapse and bounce 

When evolving the initial models described in the pre- 
vious section, the core starts to collapse. Fig. Q] shows 
the evolution of the central density for model 5B01 as a 
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FIG. 1: Evolution of the central density for the collapse model 
2501 using a semilogarithmic scaling. During the collapse the 
central density increases by 4.5 orders of magnitude. When 
reaching supra-nuclear densities, the collapse is stopped as a 
consequence of the stiffening in the hybrid EoS at about 40 
ms after the collapse was initiated. The central density finally 
approaches a new equilibrium supra-nuclear value. Shortly 
after bounce, oscillations appear in the central density (see 
inset). 



function of the Bondi time ub- The lapse of Bondi time 
as seen by an observer at infinity is defined by 



duf 



2H 



(39) 



where H = lim^-^cxD P- The conformal factor uj relates 
the two-geometry of the Bondi metric 



,^^d9^ 



to the two-geometry of a unit sphere 



ds% 



del 



sm ubc 



(40) 



(41) 



as dsl — 



u) ds . When the central density reaches nu- 
clear density at a Bondi time of about 40 ms, the pres- 
sure increases strongly according to Eq. H21(l . The central 
density grows further, but its increase is soon stopped. 
Afterwards, it drops below its maximum value, finally ap- 
proaching a quasi-equilibrium supra-nuclear value when 
a "proto-neutron star" has formed in the central re- 
gion ^26]. 

Fig. 121 shows a spacetime diagram for the core collapse 
simulation of model 2t (the main aspects are similar for 
all our models). The diagram shows different mass shells 
and the location of the shock front (thick solid line). In 
order to localize the shock front, we search for coordi- 
nate locations where the x-component of the 4-velocity 
fulfills uf — uf_^_l > s, with s being a threshold value 
for a velocity jump to be adapted (typical values for our 
simulations are s = 10-^..10-*). In addition, to com- 
pute the mass inside a fixed radius, we make use of the 
relation 



M = 47r 



-2/3 



T dr 



(42) 
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FIG. 2: Spacetime diagram for the collapse model 21. Plotted 
is the lapse of proper time as a function of the radial coordi- 
nate r. The black solid lines correspond to a subset of the light 
curves by which we foliate the spacetime (there is one light 
cone after every 5 ms, where time is measured by an observer 
at the origin). The dashed curves correspond to different mass 
shells: M = O.2M0, 0.4M0 , O.GMq, O.SM©, I.OMq, 1.2Mq . 
After about 40 ms, a shock (thick solid line) forms in the 
interior region close to the origin. The diagram was obtained 
from a global simulation with 800 radial zones, extending the 
grid to future null infinity. 



valid for the spherical collapse model 21. Fig. |21 shows 
that at the beginning of the collapse phase, the space- 
time metric is close to the Minkowski metric, which is 
reflected in the diagram by the light cones being almost 
parallel straight lines. The effects of curvature can be 
most strongly seen close to the origin (r = 0) after about 
40 ms, when the proto-neutron star has formed. We 
observe a redshift factor e^^ relating the lapse of local 
proper time at the origin to the lapse of proper time at 
infinity of ~ 1.12. 

Correspondingly, Fig. |31 shows different snapshots of 
the radial velocity at evolution times close to bounce. 
In the inner region (the so-called homologous inner core) , 
the infall velocity measured as function of radius is pro- 
portional to the radius. The homologous inner core 
shrinks with time. The outer limit of the homologous re- 
gion, i.e. the sonic point, where the local sound speed has 
the same magnitude as the infall velocity, finally reaches 
a radius of less than 10 km after about 40 ms. At that 
time, the shock front forms, which moves outwards with 
a speed of ~ 0.1c initially. During its propagation it is 
gradually slowed down by the interaction with the in- 
falling material in the outer region. It is worth to stress 
the ability of the code to resolve the steep shock front 
within only a few grid zones (typically three). This can 
be further seen in Fig.01 where we plot the rest mass den- 
sity p at the shock front for a simulation of the collapse 
model £01. 

Matter falling through the outward propagating shock 
is heated substantially. This can be seen in Fig. |31 where 
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FIG. 3: Snapshots of radial velocity profiles , plotted as 
function of radius r for the collapse model 21. The snapshots 
are taken between m_b = 30 ms and us = 45 ms, with a 
delay of 1 ms between subsequent outputs (the solid lines 
correspond to ub ~ 30,35,40,45 ms). The shock formation 
takes place at about 40 ms. In the outer part of the plotted 
region, the infall velocity of matter increases monotonically 
with time. 




FIG. 4: Surface plot of the rest mass density distribution p 
around the shock front for the collapse model £01. 50 ms after 
the collapse was initiated, the shock has reached a radius of 
about 250 km. We plot every radial zone using a radial grid 
r = lOOtan(^x) with 450 radial zones. The shock front is 
resolved with only three radial zones. The aspherical nature 
of the data is most prominent at the shock front. 



we plot the internal energy distribution e in the central 
region shortly after bounce. The figure further shows 
the contribution to the internal energy from the poly- 
tropic part, Eq. J^^J, and the thermal part, Eq. 
In the very central region, the polytropic contribution 
constitutes the dominant part. In contrast, the thermal 
energy dominates the total internal energy in the post- 
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FIG. 5: Radial distribution of the internal energy e (solid line) 
shortly after bounce {ub = 41 ms) for the collapse model 
21. The different contributions from the polytropic part ep 
(dashed line) and the thermal part eth (long-dashed line) to 
the total internal energy are also shown. In front of the shock 
which is located at a radius of ~ 45 km, the thermal energy 
vanishes. 

shock region for radii larger than a certain value (the 
shock forms off center), ~ 13 km in the specific situation 
shown in Fig. |S1 We have verified that the global energy 
balance (see Ref. j3 for rnore details) is well preserved 
in our simulations (maximum errors are of the order of 
0.5 - 1%). 

Fig. El shows two-dimensional contour plots illustrat- 
ing the dynamics during collapse and bounce for model 
*801. For this particular simulation we used a resolution 
{Nx,Ny) = (600,12). The figure displays isocontours of 
the rest mass density covering only the inner part of the 
iron core up to a radius of 30 km at 40 ms (i.e. at bounce; 
top panel), at 45 ms (when the shock has reached a ra- 
dius of 140 km; middle panel) and at 50 ms (when 
the shock wave is located at r ~ 250 km; bottom panel). 
The velocity vectors overlayed onto the contour plots are 
normalized to the maximum velocity in the displayed re- 
gion. During the collapse phase until bounce at nuclear 
densities (upper panel), the initial aspherical contribu- 
tions do not play a major role - the radial infall velocities 
dominate the dynamics. After bounce (middle and lower 
panel) the newly formed neutron star in the central re- 
gion shows nonspherical oscillations, with fluid velocities 
up to about 2 X lO'^^c. Qualitatively, the dynamics for 
the collapse model £01 is very similar to what is shown 
in Fig.Elfor model 0501. However, the particular form of 
the non-spherical pulsations created after bounce differs. 



B. Fluid oscillations in the outer core 

When analyzing the dynamical behavior of the fluid af- 
ter bounce, we find that the meridional velocity oscillates 
strongly in the entire pre-shock region. This can be seen 
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FIG. 6: Contour plot of the rest mass density distribution 
for model ©01 at a Bondi time us = 40 ms (upper panel), 
u_B = 45 ms (middle panel) and us = 50 ms (lower panel), 
obtained from a global evolution extending the grid to future 
null infinity. We only show a fraction of the core up to a 
radius of 30 km. Overlayed are velocity vectors. At bounce 
(upper panel), the matter distribution is, to a great extent, 
spherically symmetric . In the later phases (middle and lower 
panels) , the fluid dynamics is characterized by aspherical flows 
related to the oscillations of the newborn neutron star. The 
matter flow shows reflection symmetry with respect to the 
equator, which is inherent to the initial data and well pre- 
served during the evolution. 



from the solid curve of Fig. where we plot the merid- 
ional velocity component vi = rvP for model SOI as a 
function of the Bondi time, and at coordinate location 
r = 833 km and y = 0.5. These oscillations are created 
directly after the formation of the proto-neutron star in 
the central region of the numerical domain. The only pos- 
sibility to propagate information instantaneously (i.e. on 
a slice with constant retarded time u) from the central 
region to the outer layers of the iron core is through the 
metric, since sound waves would need several 10 ms to 
cover the distance. There are two possible explanations 
for these oscillations. Either they are created when gravi- 
tational wave energy is absorbed well ahead of the shock, 
or they are created by our choice of coordinates, i.e. they 
are gauge effects. In the latter case, the oscillations would 
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FIG. 7: Meridional velocity component as a function of Bondi 
time at the flxed location r = 833 km and y — 0.5 for model 
2501. The radial location was chosen well ahead of the shock. 
The solid line corresponds to the meridional velocity as ex- 
tracted in our coordinate system, V2 — ru^ , in units of the 
speed of light c. The dashed line corresponds to the merid- 
ional velocity evaluated in inertial Bondi coordinates deflned 
at future null inflnity. See text for more details. 



not be caused by a real flow, but as a consequence of the 
underlying coordinate system in which we describe the 
flow. 

To clarify the origin of the oscillations we estimate in 
the following the kinetic energy of the oscillations, as- 
suming that they are a physical effect. The average am- 
plitude of the oscillation is of the order oi V2 = 2 x 10~'*c. 
Note that V2 vanishes at the polar axis and at the equa- 
tor, so that the average velocity is substantially smaller 
than that shown in Fig. |7| Taking into account that 
the total mass in the pre-shock region is of the order 
of Mps ^ ^Mq, the kinetic energy of the oscillations is 
roughly 



kin 



^Mp,(«2)2 



2 X IQ-'^Mp.c^ 



(43) 



This energy is comparable to the total energy radiated in 
gravitational waves in a typical core collapse event 0, 0] ■ 
Transferring such an amount of energy to the pre-shock 
region seems unphysical, as gravitational waves interact 
with matter only very weakly. Instead, as we describe 
next, we conclude that the oscillations are mainly intro- 
duced by our choice of coordinates. 

Following the work of Bishop et al. 23 inertial coordi- 
nates can be established at future null infinity J7^. The 
angular inertial coordinate 9b can be constructed solving 
the partial differential equation 

(du + Ude)eB = 0, (44) 

with initial data 9b{u = 0) = 8{ii = 0). Instead of 
solving Eq. 1441) directly, we determine its characteristic 
curves, 

dO 
du 



U{9,u) 



(45) 
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(u = 0) 



(46) 



along which is constant. With suitable interpolations, 
Ob can then be determined for arbitrary angles 9. 

Making use of Eq. 144|1 . it is possible to define an "in- 
ertial" meridional 4-velocity component 



dee 

de 



J+ 



(47) 



The dashed line in Fig.[7|shows the corrected ("inertial") 
meridional velocity ru^^. Remarkably, the oscillations 
have almost disappeared, which clearly shows that gauge 
effects can play a major role for the collapse dynamics in 
the pre-shock region. 



V. GRAVITATIONAL WAVES 
A. Quadrupole gravitational waves 

The common approach to the description of gravita- 
tional waves for a fluid system relies on the quadrupole 
formula p^ . The standard quadrupole formula is valid 
for weak sources of gravitational waves under the assump- 
tions of slow motion and wave lengths of the emitted 
gravitational waves smaller than the typical extension of 
the source. The requirement that the sources of gravita- 
tional waves are weak includes the requirement that the 
gravitational forces inside the source can be neglected. 
This first approximation can be extended based on Post- 
Newtonian ex pan sions (for a detailed description see the 
recent review j28j and references therein). 

In a series of papers 1^ 113 j Winicour estab- 

lished that the quadrupole radiation formula can be de- 
rived in the Newtonian limit of the characteristic field 
equations. Let Q be the quadrupole moment transverse 
to the (9, 6) direction 



Q = gV(-) (-) 

\ r J .A \ r J .B 



where 



'■/3)d^x 



(48) 



(49) 



is the quadrupole tensor and qa, A = 2,3, is the complex 
dyad for the unit sphere metric 



d6'^ + sin 



2 /) JJ,2 



'^1{AlB)dx dx 



(50) 



As usual we use parentheses to denote the symmetric 
part. For our axisymmetric setup, Eq. (|48|l reduces to 



Q ^ TTsin^ 9 ( dr' ( sin 9' d9' r"^ pC^ cos^ 9' -]- 
Jq Jo ^2 2 

_(51) 

On the level of the quadrupole approximation the 
quadrupole news Nq reads 



Nn 



du% 



(52) 



With our null foliation it is natural to evaluate the 
quadrupole moment l|51|l as a function of retarded time, 
i.e., for the evaluation of the integral we completely relax 
the assumption of slow motion. 

It is well known 33] that the third numerical time 
derivative appearing in Eq. H52(l can lead to severe nu- 
merical problems resulting in numerical noise which dom- 
inates the quadrupole signal. Therefore, we make use of 
the fluid equations in the Newtonian limit to eliminate 
one time derivative. Deflning the "Newtonian velocities" 



vi = w'' = ^u'^, 
dx 



V2 — ru 



smf 



(53) 
(54) 



the quadrupole radiation formula 1)52(1 can be rewritten 
with the use of the continuity equation as the so-called 
first moment of momentum formula 



dU^g 



Us\n^9 [ dr' [ sm9'd9'r'^ 
^ Jo Jo 



p{vi{3cos^9' - 1) - 3t)2sin6''cos6'')). (55) 



No 



We henceforth work with Eqs. (|52|l and (|55|l for esti- 
mating the quadrupole radiation. In addition, following 
earlier work 0, we deflne the quantity A2Q, which 
enters the total power radiated in gravitational waves in 
the quadrupole approximation as 



dE 

dUB 



1 fdA^i\ 
327r \ duB ' 



(56) 



A2Q also arises as coefficient for the quadrupolar term in 
the expansion of the quadrupole strain (i.e. t he g ravita- 
tional wave signal) /i+ in spherical harmonics j34| 



h+{uB) 



1 /15 



R ' 



(57) 



where R denotes the distance between the observer and 
the source, 
moment as 



A20 can be deduced from the quadrupole 



\E2 

^20 



16 3 d^ 
7r2 



du^ 



dr' 



sin 9'd9'r 



cos2 9'-- 



(58) 



or alternatively using the first moment of momentum for- 
mula in order to eliminate one time derivate, in analogy 
to the transition from Eq. 1)52(1 to Eq. ((55(1 . i.e., 



aE2 

^20 



16 



7r2 ■ 



R 

dr' 
^0 



sin 6''d6'V 



15 duB 

X p(i;i(3cos2 6'' - 1) - 3t;2sine''cos6l') 



(59) 



As shown in Fig.|Slwe find good agreement when com- 
puting the wave strain A2Q using Eq. ((58(1 and Eq. 1(59(1 . 
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FIG. 

the collapse model 2301. The solid curve shows the result 
using the first moment of momentum approach Eq. I59II . the 
dashed line is based on Eq. 1581 . The good agreement found 
between both approaches shows that our general relativistic 
fluid evolution is internally consistent. 



In order not to have the time derivatives dominated by 
numerical noise, we have averaged the matter contribu- 
tion in the integrands of Eq. H58|) and Eq. (|59f) over a 
few neighboring grid points before calculating the time 
derivatives. 

This result checks the implementation of the conti- 
nuity equation and, as this equation is not calculated 
separately but as a part of a system of balance laws, 
it also checks the overall implementation of the fluid 
equations in the code. We note that the equivalence 
between Eq. H58|) and Eq. H59|l is only strictly valid in 
the Minkowskian limit and for small velocities, which is 
the origin for the observed small differences between the 
curves in Fig. |S1 Substituting p by pu^e^^ in Eq. 
and by pe^^ in Eq. (|59|l . by which we restore the equiva- 
lence in a general relativistic spacetime, we find excellent 
agreement between the two approaches for calculating 

AE2 
^20 ■ 

Since we are imposing only small perturbations from 
spherical symmetry, we expect a linear dependence of 
the non-spherical dynamics and the gravitational wave 
signal as a function of the perturbation amplitude. We 
have verified in a series of runs that the amplitude of the 
quadrupole moment (and thus the quadrupole radiation 
signal) indeed scales linearly with the amplitude of the 
initial perturbations (see Fig.|5Jl. This observation marks 
another important test for the correctness of the global 
dynamics of our code. 

On the other hand, when comparing the quadrupole 
news defined in Eq. H52|l or Eq. (|55|l with the Bondi news 
signal N evaluated at future null infinity (which is defined 
in Eq. H68() below), we find important discrepancies. This 
can be seen in Fig. ^| where we plot both, the Bondi 
news and the quadrupole news for model 5801. We note 
that the differences manifest themselves not only in the 



FIG. 9: Quadrupole moment Q (in units c = G = Mq = 
1) as a function of time for three models of type 53 with 
perturbation amplitude A = 0.01, A = 0.05 and A — 0.1. 
The first two results are rescaled with respect to A = 0.1 
assuming a linear dependence. All three curves overlap in the 
diagram. The quadrupole moment (and hence the quadrupole 
signal) scales linearly with the amplitude of the perturbation 
in the chosen parameter region. 
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FIG. 10: Bondi and quadrupole news as a function of time for 
model 2301. The solid curve corresponds to the quadrupole 
news according to Eq. H55|l . the dashed curve to the Bondi 
news signal. For visualization reasons, we have divided the 
quadrupole news result by 50. Remarkable disagreement is 
found between both signals. 



amplitude of the oscillations, but also in the frequencies 
of the signals. This behavior is clearly different from the 
one we observed in the studies of neutron star pulsation 
carried out in Ref. i where both signals showed very 
good agreement. 

As mentioned above, the quadrupole formula is only 
the first term in a Post-Newtonian expansion for the grav- 
itational radiation. The next, non- vanishing contribution 
to the gravitational strain for our axisymmetric configu- 
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ration is the hexadecapole contribution, which reads 



800 



/^r-^/f-n^^a-Zsin^^)^. (60) 



The quantity Afg is defined as 



aE2 _ " A,fE2 



(61) 
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X p(7cos*6i'-6cos^6l' + -), (62) 
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or alternatively 



d3 



aE2 _ tyE2 

iVf 2 ^ l}/5^i Z''^^/ r sine' dO'r'' 
63 7o Jo 



(63) 



X p(i;i(7cos^6l'-6cos2 6i' + -) 

+W2(3-7cos2 6'')sin6''cos6''). (64) 



By extracting the hexadecapole moment for the above 
result, we found, however, that the associated amplitude 
is too small in order to explain the observed differences 
in Fig. 1101 In addition, one would expect in general that 
the contribution of the hexadecapole moment increases 
the amplitude of the approximate signal. However, the 
amplitude of the quadrupole news in Fig. 1101 is already 
much larger than that of the Bondi news evaluated at 

As we discussed in the preceding section, the global 
dynamics of the core collapse and bounce is correctly 
reproduced with our numerical code (see also the vali- 
dation tests in the Appendix). We have strong evidence 
that the quadrupole signals extracted from our collapse 
simulations do not correspond to physical gravitational 
wave signals. In the following, we describe the different 
arguments which support this claim. 

First, if the quadrupole radiation signal corresponded 
to the true physical signal, it would be very difficult 
to understand why the Bondi signal has a significantly 
smaller amplitude. In the calculation of the Bondi news, 
Eq. H68|l . the contribution of the different terms are rel- 
atively large and add up to a small signal (see below). 
Under the assumption that the quadrupole news signal 
is correct and the Bondi news signal is wrong, it is ex- 
tremely unlikely that possible errors in the contribution 
to the Bondi news add up to a very small signal. 

Second, we have performed comparisons between our 
numerical code and the code of Refs. finding much 

larger amplitudes for the quadrupole gravitational wave 
signal in our case. However, we note that comparing 
the results of both codes in axisymmetry is ambiguous, 
as possible differences might have different explanations. 
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FIG. 11: Radial contribution to the quadrupole mo- 
ment. We plot the value of the integral Q{r) = 
TT sin'^ e dr' sine' de'r'* cos'^ e' - i) as a function of 
the radial coordinate r for different values of time. The data 
is plotted after a fixed number of time steps, starting with 
initial data at us = ms (upper solid curve). The data was 
taken from a simulation of model 2301. Large amplitude os- 
cillations of the quadrupole moment, as they can be seen in 
Fig. El can only be created - at least shortly after bounce - 
in the outer region of the infalling matter well in front of the 
shock. 



For example, the use of the conformally flat metric ap- 
proach in y| is clearly an approximation to general 
relativity, which should create some differences. Fur- 
thermore, the coordinate systems used in both codes for 
the computation of the quadrupole moment are different. 
Only in our code, the quadrupole moment is evaluated 
on a light cone, i.e. as a function of retarded time. 

A third and physically motivated argument stems from 
the spatial distribution of matter in our simulation. As 
it can be seen from Fig. 1111 the main contribution to 
the radial integral of the quadrupole moment comes from 
the outer, infalling layers of matter. These outer layers 
are responsible for the oscillations in the quadrupole mo- 
ment, which can be seen in Fig. |^ Following the same 
reasoning as in the previous section it is obvious to con- 
clude that the calculation of the quadrupole moment is 
also affected by our choice of coordinates, i.e. by gauge 
effects. 

For all these reasons we extract the quadrupole mo- 
ment in the angular coordinate system defined by 
Eq. (|44|l . However, introducing the inertial angular co- 
ordinate does not help to obtain a better agreement 
between quadrupole and Bondi signals, the extracted 
quadrupole moment almost agrees with the results shown 
in Fig.ini Since the difference of Bondi time between the 
different angular directions on our Tamburino-Winicour 
foliation is in general of the same order as the lapse of 
time for one time step, we expect a similar result when 
evaluating the quadrupole moment at a fixed inertial 
time. However, by prescribing the necessary coordinate 
transformations to define Bondi coordinates only at J^~^ , 
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we do not take into account an inertial radial coordinate, 
which should be used for the evaluation of the quadrupole 
moment. 

As already mentioned before, in Ref. jT^ we found 
good agreement between the Bondi signal and the 
quadrupole signal when calculating gravitational waves 
from pulsating relativistic stars. Hence, the obvious ques- 
tion arises of why the quadrupole formula could be ap- 
plied in those scenarios. The answer lies in the small 
velocities encountered in the problem of neutron star pul- 
sations. Whereas the typical maximum fluid velocities in 



the oscillation problem are of the order of 10 



10- 



fluid velocities of up to 0.2c are reached for the core col- 
lapse scenario. Furthermore, due to the non-spherical 
dynamics of the proto-neutron star formed in the inte- 
rior of the collapsed region, the metric can pick up gauge 
contributions which are created as a consequence of our 
requirement to prescribe a local Minkowski frame at the 
vertex of the light cones. Gauge contributions may also 
play a more important role in the collapse scenario due to 
the enlarged radial extension of the collapsing iron core 
(about 1500 km), which is much larger than the corre- 
sponding one for neutron star pulsations (about 15 km). 

We note that since the collapse involves fluid velocities 
of up to 0.2c, it is not obvious whether the functional 
form for the quadrupole moment established in the slow 
motion limit on the light cone will still be valid. In fact, 
the situation could be similar to the case of the total mass 
of spacetime, where a naive definition, even in spherical 
symmetry, as 



An 



(65) 



would only be a valid approximation for small fluid ve- 
locities. This can be understood from the comparison 
with the expression of the Bondi mass in the form 



Mf 



All / r'2[p(l-|-e)(-u-7i„)-p(l+u«u„)]dr', (66) 
Jo 



(no summation is involved in this expression). Only van- 
ishing fluid velocities, i.e. u^Uu = —1, ensure that the 
two masses are equal, Mn = Mb- 

We experimented with possible alternative functional 
forms for the quadrupole moment which result in signif- 
icant differences. An unambiguous clarification of which 
functional form has to be used for the quadrupole mo- 
ment in the extended regime of validity of large fluid 
velocities could only be obtained by a derivation of the 
quadrupole formula in the Tamburino gauge. However, 
technical complications for such a derivation are so se- 
vere that it has only been accomplished for a simplified 
radiating dust model js^f (see the related discussion in 
Ref. 111). 



B. The Bondi news signal 

The numerical extraction of the Bondi news is a very 
complicated undertaking. Reasons for possible numerical 
problems are diverse: First, its extraction involves calcu- 
lating non-leading terms from the metric expansion at 
future null infinity. All the metric quantities are global 
quantities, and are thus sensitive to any numerical prob- 
lem in the entire computational domain. Second, when 
calculating the gravitational signal in the Tamburino- 
Winicour approach, one has to take into account gauge 
effects. For the present calculations of the gravitational 
wave signal from core collapse, the gauge contributions 
are indeed the dominant contribution, which can easily 
influence the physical signal. 

We have described in detail the formalism and numer- 
ical methods to deal with gravitational waves without 
approximation in our axisymmetric characteristic code 
in Ref. In the following, we will only repeat the 

most important aspects. The total energy emitted by 
gravitational waves to infinity during the time interval 
[u, u -\- du] in the angular direction [y, y -f dy] is given by 
the expression 



dE = ^N^uj^e^^dydu, 
where the Bondi news function N reads 



(67) 
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(68) 



K, c, H and L are defined by a power series expansion of 
the metric quantities at J'^ as follows. 



7 = K+-+0{r-^), 
r 

U = L + 0{r-^). 



(69) 

(70) 
(71) 



We plot in Fig. ^| the different contributions to the 
Bondi news for the collapse model 5B01. It becomes clear 
from this plot that a very accurate determination of the 
metric fields is essential. As it can be further seen in this 
figure, the metric quantities show high frequency numer- 
ical noise, as soon as the shock forms (at a Bondi time of 
about 40 ms). In order to demonstrate that the noise is 
actually created at the shock, we plot in Fig. ^| the lo- 
cation of the shock together with the gravitational wave 
signal. Clearly, the noise is created by the motion of the 
shock across the grid, its temporal behavior following the 
discontinuous jumps of the shock between adjacent grid 
cells. We note that due to the coarser radial resolution 
used in the outer layers of the core, the frequency of the 
noise slowly decreases with time. 

As we have pointed out in the previous section, the 
shock front is well captured in only a few radial zones 
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FIG. 12: DiflFerent contributions to the Bondi news. The solid 
curve corresponds to the term involving Cu (first addend) in 
Eq. 1)68^ . the dashed curve to the contribution from the sec- 
ond and third addend. By summing up both contributions 
we obtain the Bondi news, which is close to zero. In addition, 
we note that when separating the third addend into angu- 
lar derivatives of H and ui, each single contribution has an 
amplitude 23 times larger than what is shown in the figure. 



with our high-resolution shock-capturing scheme. It 
might seem surprising that a small localized error cre- 
ated in a few radial zones can have such a large effect on 
the gravitational wave signal. However, one has to keep 
in mind that the radial integration of the metric vari- 
ables picks up this error and propagates it to future null 
infinity instantaneously. It is important to stress that 
the effect of the numerical noise on the dynamics of the 
collapse and bounce is entirely negligible. However, the 
extraction of the Bondi news signal is extremely sensitive 
to it. 

We have verified that the frequency of the noise in- 
creases, as expected, with radial resolution. Unfortu- 
nately, its amplitude does not decrease substantially with 
radial resolution, at least not in the resolution regime 
accessible to us [IHl- Therefore, we tried to eliminate 
the noise by different methods. In a first attempt, we 
smoothed out the shock front, either in the hydrodynam- 
ical evolution itself or before using the fluid variables in 
the source terms of the metric equations. In both cases, it 
was impossible to obtain a smooth signal without chang- 
ing the dynamics. In a second attempt, following the 
work of |37l |. we rearranged the metric equations elimi- 
nating second derivatives which might be ill-behaved at 
the shock. Defining a metric quantity 



(72) 



and solving the hypersurface equations successively for (3, 
X, U and S, it is possible to eliminate all second deriva- 
tives from the hypersurface equations. Unfortunately, 
the noise is not significantly reduced by this rearrange- 
ment of the metric equations. Finally, going to larger 
time steps for the fluid evolution only - solving the met- 
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FIG. 13: Upper panel: Bondi news as a function of time. 
High frequency noise is overlayed on top of a small frequency 
modulation. Lower panel: Time evolution of the radial loca- 
tion of the cross section of the shock front with the equator. 
Due to the finite resolution, the location of the shock wave 
moves discontinuously. The frequency of these jumps coin- 
cides with that of the noise in the Bondi news. Once created 
at the shock, the noise is propagated instantaneously to in- 
finity through the numerical solution of the metric equations. 



ric equations several times between one fluid time step - 
was not effective either. 

After these attempts we decided to eliminate the noise 
from the gravitational wave signals only after the nu- 
merical evolution. We experimented with two different 
smoothing methods. In the first method, we calculate 
the Fourier transform of the data, and eliminate all fre- 
quencies beyond a certain threshold frequency (of about 
5-10 kHz). Then, when transforming back from Fourier 
space all the high-frequency part of the data is removed. 
In a second method we simply average the signal over 
a few neighboring points. We have applied this second 
method in what is described below. 

Fig. ^] shows the Bondi news signal for the collapse 
model QSOl. The figure focuses on the part of the signal 
around bounce. After the initial gravitational wave con- 
tent is radiated away (in the first 5 ms, not depicted in 
the figure), the signal in the collapse stage is very weak. 
This is expected, as the dynamics is well reproduced by 
a spherical collapse model during this stage. At bounce, 
the Bondi news shows a spike. Afterwards, a complicated 
series of oscillations is created due to the pulsations of the 
forming neutron star and the outward propagation of the 
shock. Typical oscillation frequencies are of the order of 
0.5— 1 kHz, at which the current gravitational wave laser 
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FIG. 14: Bondi news as a function of Bondi time for the col- 
lapse model 2301. The displayed time interval covers the late 
collapse stage until several ms after bounce at about t = 40 
ms. During the collapse stage, the gravitational wave signal 
is negligible. After bounce a complicated series of oscillations 
sets in. 
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FIG. 15: Bondi news as a function of time for the collapse 
model £01. The bounce at about 41 ms is characterized by a 
large spike in the gravitational wave signal. After bounce, the 
signal shows oscillations, with a principal frequency of about 
0.35 kHz. 



interferometers have maximum sensitivity. 

Correspondingly, Fig. 1151 shows the Bondi news signal 
for the collapse model £01. Here again, after radiating 
away the initial gravitational wave content, the collapse 
phase is characterized by very small radiation of gravita- 
tional waves. At bounce, we again observe a strong spike 
in the signal. Afterwards, the oscillations in the signal 
are rather rapidly damped. 

We stress that as a consequence of the necessary 
smoothing techniques applied, only the main features of 
the gravitational wave signals in Figs. 1141 and 1151 are re- 
liably reproduced. This also applies to possible offsets of 



the Bondi news, which affect in particular the gravita- 
tional wave strain. Comparing the Bondi news function 
for the different collapse models of type *B, we observe 
to good approximation a linear dependence of the Bondi 
news with the perturbation amplitude. This is reflected 
in the total energy radiated away in gravitational waves, 
which scales quadratically with the amplitude of the ini- 
tial perturbation. A summary of the results on the grav- 
itational wave energy is listed in in Tabled 



TABLE I: Total energy radiated in gravitational waves during 
the first 50 ms for the collapse simulations of type 23. The 
initial gravitational wave content is the dominant contribution 
to the total energy. This energy scales quadratically with the 
amplitude of the initial perturbation, as can be inferred from 
the last column, where the corresponding energies have been 
rescaled with respect to that of collapse model 2301. 



model 


total energy radiated [Mq] 


rescaled result [Mq] 
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4.31 X 10"^ 
1.08 X 10"^ 

4.32 X lO"'^ 


4.31 X 10"^ 

4.32 X 10"^ 
4.32 X 10"^ 



VI. DISCUSSION 

We have presented first results from axisymmetric core 
collapse simulations in general relativity. Contrary to 
traditional approaches, our framework uses a foliation 
based on a family of light cones, emanating from a reg- 
ular center, and terminating at future null infinity. To 
the best of our knowledge, the characteristic formulation 
of general relativity has never been used before in simu- 
lations of supernova core collapse and in the extraction 
of the associated exact gravitational waves. Our axisym- 
metric hydrodynamics code is accurate enough to allow 
for a detailed analysis of the global dynamics of core col- 
lapse in general. But we have not found a robust method 
for the (Bondi news) gravitational wave extraction in the 
presence of strong shock waves. 

Comparing our results to other recent work on rela- 
tivistic supernova core collapse 0, ) it is not surprising 
that numerical noise in the gravitational waveforms is 
more noticeable in our approach. Whereas in the con- 
formal flat metric approach employed in 0, Q the met- 
ric equations of general relativity reduce to elliptic equa- 
tions, which naturally smooth out high-frequency numer- 
ical noise, we solve for the gravitational wave degrees of 
freedom directly using the full set of field equations of 
general relativity, and hence we have to solve a hyperbolic 
equation. It remains to be seen whether a similar numer- 
ical noise to the one we find when extracting the gravi- 
tational wave signal will be encountered in core collapse 
simulations solving the full set of Einstein equations in 
the Cauchy approach. In this respect we mention recent 
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axisymmetric simulations by Shibata using a conformal- 
traceless reformulation of the ADM system ||38|j where, 
despite of the fact that long-term rotational collapse 
simulations could be accurately performed, gravitational 
waves could not be extracted from the raw numerical 
data since their amplitude is much smaller than that of 
other components contained in the metric and/or numer- 
ical noise. 

With the current analysis we have presented in this 
paper, it is not obvious how the numerical noise of the 
Bondi news can be effectively eliminated. Including rota- 
tion in the simulations, which would be the natural next 
step for a more realistic description of the scenario, could 
help in this respect. Due to the global asphericities intro- 
duced by rotation, one would expect, in general, gravita- 
tional wave signals of larger amplitude, which could make 
the numerical noise less important, if not completely ir- 
relevant. In addition to this possibility we propose the 
following methods to improve the extraction of the grav- 
itational wave signals: In a first approach one should try 
to rearrange the metric equations by introducing auxil- 
iary fields which could effectively help to diminish the im- 
portance of high-order derivatives, especially of the fluid 
variables, which can be discontinuous. Unfortunately, 
to the best of our knowledge, there is no clear guide- 
line to what is really needed to eliminate the numerical 
noise completely, apart from the hints given by [37j . Our 
attempts in this direction have not yet been successful, 
but we believe there is still room for improvement. Al- 
ternatively, one should try to implement pseudospectral 
methods for the metric update. Pseudospectral methods 
would allow for a more efficient and accurate numerical 
solution of the metric equations. In a third promising 
line of research we propose to consider the inclusion of 
adaptive grids and methods of shock fitting into the cur- 
rent code. With the help of an adaptive grid, one could 
try to arrange the entire core collapse simulation in such 
a way that the shock front always stays at a fixed lo- 
cation of the numerical grid. By avoiding the motion 
of the shock front across the grid, one would expect the 
noise in the gravitational wave signals to disappear. But 
already increasing the radial resolution substantially in 
the neighborhood of the shock front could help to obtain 
an improved representation of the shock. All these issues 
are ripe for upcoming investigations. 
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APPENDIX A 

In this Appendix we present tests specifically aimed 
to calibrate our code in core collapse simulations. The 
reader is addressed to Ref. 18] for information on fur- 
ther tests the code has successfully passed concerning, 
among others, long-term evolutions of relativistic stars 
and mode-frequency calculations of pulsating relativistic 
stars. 



1. Shock reflection test 

In order to assess the shock-capturing properties of 
the code, we have performed a shock reflection test in 
Minkowski spacetime. This is a standard problem to cal- 
ibrate hydrodynamical codes j.39j . A cold, relativistically 
inflowing ideal gas is reflected at the origin of the coor- 
dinate system, which causes the formation of a strong 
shock. We start the simulation with a constant density 
region, where p = po, u"^ = u^^ and e = er = (we 
set e « 10~^^ for numerical reasons). From the continu- 
ity equation it follows that the rest mass density in the 
unshocked region obeys 



PR{u,r) = po 1 



(Al) 



From momentum conservation arguments, it is clear that 
the velocity in the shocked region vanishes, w£ — 0. Eval- 
uating the Rankine-Hugoniot jump conditions for the 
fluid equations, we obtain: 
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Here, s denotes the shock speed and ps = Pr{u, r 
the rest mass density in front of the shock. 

We performed this test with different values of the 
fluid velocity, and different schemes for the fluid evolu- 
tion. Fig. ^] shows the results for an ultrarelativistic 
flow {u^ = —0.9999c). For this particular test we used 
the HLL solver and increased the numerical viscosity by 
a factor 2 in order to damp small post-shock oscillations. 
The agreement with the analytic solution is satisfactory, 
and the shock front is very steep, being resolved with 
only one or two radial zones. The deviation close to the 
origin is a well-known failure of flnite-difference schemes 
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FIG. 16: Shock reflection test for an ultrarelativistic flow with 
vT = -0.9999c and po = 8 and EoS p = 5 x 10~*pt , which 
is reflected at the origin of the coordinate system. We have 
plotted different fluid quantities at an evolution time u = 
2.029 as a function of the radial coordinate r. Top panel: 
fluid velocity . Middle Panel: pressure p. Bottom Panel: 
rest mass density p. The solid line corresponds to the exact 
solution, the crosses are taken from our numerical simulation. 
For the above result, we made use of a non- equidistant radial 
grid r = x / (1 — x^) with 800 radial zones, the MC slope 
limiter and the HLL approximate Riemann solver. 

for this problem (see, e.g. |40|). which is not important 
for our purposes. 



2. Convergence tests 

We describe now some tests which check various prop- 
erties of spherically-symmetric core collapse. We choose 
a particular collapse model, for which the initial central 
density is pc — 1-62 x 10^* (in units G = c = Mq = 1), 
the polytropic constant is k = 0.46, and the collapse is 
induced by resetting the adiabatic exponent to Fi — 1.3 
(for the equilibrium model with F — |). We use the 
hybrid EoS discussed in Section III CI 



a. Thermal energy during the infall phase 
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FIG. 17: Thermal energy as a function of the radial coordi- 
nate X at 30 ms for a compactified grid r = for different 
resolutions. Due to numerical errors, the thermal energy is 
different from zero. Deviations converge to zero, the conver- 
gence rate is 2. 
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FIG. 18: Evolution of the central density for a core collapse 
induced by resetting the adiabatic exponent to Fi — 1.30. 
The central density increases by almost 5 orders of magnitude, 
before the core bounces. Afterwards the central density stays 
almost constant. The different lines correspond to different 
grid functions and resolutions, see run 4, 1 and 5 in Table IHl 



Before the central density of the collapsing core 
reaches nuclear densities, the collapse is exactly adia- 
batic. Hence, the thermal energy, which vanishes ini- 
tially, should vanish throughout this phase. This can 
be easily checked and used for convergence tests. Fig.lTTI 
shows the result after an integration time of 30 ms (when 
the central density has increased by roughly a factor 10). 
We find that the errors from the exact result eth — 
converge to zero, the convergence rate is 2. Note that 
although eth ^ from the physical point of view, the 
numerical errors can result in negative values for eth. 



b. Time of bounce 

Using the axisymmetric code developed by Dim- 
melmeier et al. 0! 0, based on the conformally flat 
metric approach, we can perform comparisons between 
the evolutions of the same initial models. As the con- 
formally flat metric approximation is exact for spherical 
models, comparisons in spherical symmetry are unam- 
biguous. 

We define the time of bounce as the time, when the 



central density reaches its maximum. In order to start 
with the same initial data we initiate the collapse by ray- 
tracing the evolution of Dimmelmeier's code to obtain 
the initial data on our null cone. There is no principal 
advantage in starting with initial data on a null cone or 
on a Cauchy slice. Ideally, results from stellar evolution 
would give exact initial conditions for the core collapse, 
thus eliminating the artificial procedure of resetting T 
to initiate the collapse. Fig. ^] shows the evolution of 
the central density for the relativistic code of Q and the 
results of our null code for two different grid functions. 
Table Ull summarizes our results for the time of bounce. 

Assuming our code is exactly second order convergent 
and extrapolating our results to an hypothetical infinite 
resolution, we obtain from runs 3 and 4, that the infinite 
resolution run bounces after 38.65 ms. This is internally 
consistent, a comparison of runs 2 and 4 results in a value 
of 38.66 ms. Using an even higher resolution for a differ- 
ent grid function in run 5, we observe a time of bounce 
close to the converged result. Our results on the time of 
bounce are in very good agreement with the result of , 
who find a value of 38.32 ms. The observed difference 
of less than 1% is either due to the fact that the result 
of [1| is not converged, or due to the different radial co- 
ordinates used in both codes, and thus small differences 
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in the initial data. 



TABLE II: Times of bounce for different grid functions and 
resolutions. 





code 


grid 
function 


radial 
resolution 


time of 
bounce [ms] 


1 


CFC code \S] 


see f8] 


80* 


38.32 


2 


null code 


^ ^ 150x^ 


600 


40.86 


3 


null code 


J. _ 150x^ 


800 


39.90 


4 


null code 


J. _ 150:c 


1000 


39.45 


5 


null code 


r = 100tan(f x) 


1200 


38.92 



* This number for the radial resolution cannot be directly 
compared to the values of our code, as we resolve the exterior 
vacuum region up to future null infinity with our code as well. 



As it can be seen in Fig. ^1 the comparison does not 
only give very good agreement for the time of bounce, but 
also for the dynamics of the central density in general. 
This is very important, since it shows that the global 
dynamics of the core collapse is correctly described in 
our numerical implementation. 
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